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ABSTRACT 

Waveform  multigrid  method  is  an  efficient  m<'tlio«l  lor  solving  certain  classes  of  time- 
(h'pi'udent  I’DEs.  This  paper  studi<'s  tlie  relationship  In'tween  this  method  and  the  analogous 
multigrid  metluul  for  steady-state  lU'ohh'ins.  Csing  a  Fourier- Laplace  analysis,  practical 
convergence  rate  estimates  of  the  waveform  nudtigrid  iterations  art'  ohtaiiu'd.  Experimental 
results  show  that  the  analysis  yiehls  accuratt'  performance  prt'tlict ion. 
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1.  Introduction.  Waveform  relaxation,  also  called  dynamic  iteration  or  Picard- 
Linderlof  iteration  [7],  is  a  technique  for  solving  ordinary  differential  systems  of  initial- 
value  type  [5]  [8].  Its  key  idea  is  to  solve  an  ordinary  differential  system  in  n  variables 
by  solving  sequence  of  subsystems  in  fewer  variables.  This  nature  of  the  method  allows 
independent  integration  with  different  timesteps  for  each  of  subsystems.  Thus  the 
method  is  usually  considered  in  the  context  of  parallel  or  multirate  algorithms  [2]. 

The  convergence  of  the  waveform  relaxation  methods  has  been  analyzed  by  Miekkala 
and  Nevanlinna  [7].  They  discussed  this  issue  over  infinite  time  interval  and  showed 
that,  for  linear  heat  equation  with  standard  spatial  discretization,  the  convergence  rates 
of  the  waveform  iterations  are  similar  to  those  for  the  analogous  steady-state  problems. 
Therefore  the  convergence  could  be  too  slow  for  the  waveform  relaxation  to  be  compet¬ 
itive  with  standard  timestepping  methods. 

Multigrid  techniques  (in  space)  can  be  incorporated  to  accelerate  the  convergence. 
This  was  studied  by  Lubich  and  Ostermann  [6]  who  compared  the  multigrid  performance 
of  waveform  iteration  with  liiat  of  static  iteration  (i.e.,  the  iteration  for  steady-state 
problems)  for  one-dimensional  heat  equation.  Combining  the  analysis  on  the  smoothing 
rate  of  high  frequencies,  they  conjectured  the  waveform  multigrid  performance  for  two- 
dimensional  case.  They  showed  that  the  typical  multigrid  acceleration  can  be  achieved 
with  an  estimated  convergence  rate,  which  is  similar  but  not  quite  as  good  as  the  one 
for  the  steady-state  problems. 

A  number  of  parallel  waveform  multigrid  algorithms  have  been  proposed  and  im¬ 
plemented  [9]  [11].  Analytical  and  experimental  results  have  shown  that  the  waveform 
multigrid  methods  can  be  implemented  on  a  parallel  computer  with  satisfactory  effi¬ 
ciencies. 

In  this  paper  we  study  the  relationship  between  the  waveform  multigrid  method  for 
time-dependent  PDEs  and  the  standard  multigrid  for  the  corresi)onding  steady-state 
problems.  This  study  is  important  since  the  steady-state  multigrid  has  been  investi¬ 
gated  extensively  while  the  properties  of  waveform  multigrid  algorithm  are  relatively 
unknown.  We  present  a  Fourier-Laplace  analysis  for  obtaining  practical  convergence 
rate  estimates  of  waveform  multigrid  iterations.  The  approachs  used  in  this  paper  are 
simple,  applicable  to  wide  class  of  applications,  and  provide  insight  into  the  details  of 
the  basic  interaction  between  the  coarse  grid  correction  and  the  waveform  relaxation. 

This  paper  is  organized  cis  follows.  Section  2  briefly  introduces  the  waveform  re¬ 
laxation  method.  In  Section  3,  the  waveform  relaxation  and  multigrid  iteration  are 
combined.  A  theorem  is  proved  which  indicates  that  the  convergence  rates  of  waveform 
multigrid  are  essentially  the  same  a.s  those  for  the  analogous  steady-state  problems  when 
number  of  smoothings  is  large.  Section  4  gives  the  details  of  Fourier-Laplace  analysis, 
which  is  used  for  obtaining  exact  convergence  rates.  As  an  example,  analytical  com¬ 
parison  of  multigrid  convergence  rates  in  practical  number  of  smoothings  is  made  for  a 
two-dimensional  heat  equation.  Finally,  the  comparison  to  the  measured  convergence 
rates  is  presented  in  Section  5. 

2.  Waveform  relaxation  method.  Waveform  relaxation  method  was  originally 
proposed  foi  solving  ordinary  differential  equations  consisting  of  subproblems  with  few 
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external  variables  in  circuit  simulation  [5].  Unlike  standard  time  stepping  methods, 
it  iteratively  partitions  a  large  system  into  loosely  coupled  subsystems  and  integrates 
each  subsystem  independently.  Hence,  it  is  well  suited  for  parallel  machines,  especially 
massive  parallel  machines. 

The  method  can  be  described  as  follows.  Consider  a  linear  initial  value  problem 
du 

(1)  —  +  Au  =  f,  <  >  0,  u(0)  =  uo- 

(it 

Let  A  be  split  a.s  A  =  M  —  N.  Under  certain  conditions,  Ec|.(  1 )  can  be  solved  iteratively 

by 

(2)  +  Mub')  =  +  /,  <>  0,  n<’d(o)  =  uo, 

which  is  equivalent  to  the  integral  equation 

(3)  u*")  =  +  0, 

where  «S  is  a  linear  integral  operator  on  L^(R'^ ,C’^)  ( 1  <  p  <  oo)  with  kernel  k^it)  = 
and  0  is  a  function  depending  upon  /  and  M. 

The  general  convergence  results  of  the  scheme  (2)  over  entire  interval  t  >  0  were 
given  by  [7].  In  their  work,  the  Laplace  transformation  was  applied  to  the  time  variable 
t  and  it  showed  that  the  convergence  rate  for  scheme  (2)  was  given  by  the  spectra! 
radius  of  <S,  derived  as 


(4) 


p{S)  =  meiX^p{S{z)),  S{z)  :  Laplace  transform  of  S. 


We  shall  follow  their  approach  and  concentrate  on  the  damped  .Jacobi  relaxation  and  the 
red-black  Clauss-Seidel  relaxation  for  linear  equations  with  time-independent  coefficients 
of  the  form 

du 

(5)  —  +  Lu  =  f,  <  >  0,  u(0)  =  Uo, 

where  L  is  a  linear  elliptic  operator.  Let  Lh  be  a  discrete  approximation  of  L  and 
assume  it  can  be  written  as 

(6;  Lu  =  d 


I  -B 
-R  I 


Then  the  damped  Jacobi  relaxation  and  the  red-black  Causs-Seidel  relaxation  have  the 
Laplace  transform 


(7) 


Sj{z)  =  {zI  +  Mjr'Nj  = 


(i-1)/  B 


Mj  =  -/,  Nj  =  Mj-  U, 

U) 


0  <  u;  <  1; 


•2 


(8) 


.9r;5(^’)  =  {zl  +  MGsrNcs 


Mgs  =  d 


/  0 

-R  I 


^as  —  Mgs  — 


In  the  following  sections,  the  subscripts  of  the  matrices  will  be  dropped  when  the  context 
is  clear. 

3.  Waveform  versus  steady-state  multigrid  method.  The  convergence  of  the 
waveform  relaxation  method  can  be  accelerated  if  the  multigrid  technique  is  incorpo¬ 
rated  in  space.  In  this  section,  we  shall  show  that,  the  convergence  rate  of  waveform 
multigrid  iteration  converges  to  that  of  steady-state  multigrid  iteration  as  the  number 
of  smoothings  increases.  In  next  section,  a  Fourier- Laplace  analysis  will  be  introduced 
to  show  that  the  performance  of  the  waveform  multigrid  iteration  for  time-dependent 
PDFs  is  virtually  as  good  as  the  standard  multigrid  iteration  for  the  corresponding 
steady-state  problems  for  practical  number  of  smoothings. 

The  multigrid  method  is  adapted  to  the  waveform  iteration  for  a  time-dependent 
PDE  in  the  following  w'ay.  First,  the  equation  is  discretized  in  space  to  give  a  semi¬ 
discrete  problem.  Next,  the  multigrid  iteration  with  waveform  relaxation  is  applied 
to  the  space  variables.  As  an  example,  a  two  grid  V-cycle  for  Eq.(5)  is  illustrated  as 
follows; 

•  Perform  ui  pre-smoothings: 

+  Mu^"^  =  -}-  /,  <>  0,  =  Wo, 

at 

where  Lh  =  M  —  N,  v  =  1,2,  -  is  given. 

•  Restrict  the  defect  from  grid  h  to  grid  H: 

(9)  4  :=  -  /,  dn  :=  7,^4. 

•  On  the  coarse  grid,  solve 

dv 

—  +  Lhv  =  dn,  v(0)  =  0. 

•  Correct 


where  is  a  suitable  interpolation  from  grid  H  to  grid  h. 
•  Perform  V2  post-smoothing®  on  v. 
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Similar  to  the  Full-Approximation-Srheme  (FAS)  formulation  of  multigrid,  one  can 
formulate  a  coarse  grid  problem  as 


duf^ 

dt 


"u'-’j  +  iw/ 


HJ..) 


wh(0)  =  /”  uo 

and  the  correction  step  as 

u  = 

in  order  to  handle  non-linear  problems.  This  however  should  be  used  in  an  FMG 
algorithm  where  the  problems  is  solved  first  on  coarse  levels  to  obtain  a  good  initial 
approximation  to  fine  ones. 

The  error  —  u  of  a  complete  two  grid  V'-cycle  iteration  described  above 

satisfies 

(10)  = 

for  an  integral  operator  V,  which  has  the  Laplace  transform  (see  [6]) 

(11)  C(.-)  =  .5(z)”>(/-/^(z  +  iH)-';,''(2+U))i’(2)”',  Rez>0. 

where  S{z)  is  the  Laplace  transform  of  tiie  smoother  being  used.  In  order  to  indicate 
the  dependency  of  V  and  V{z)  on  the  number  of  smoothings  v  —  U]  -I-  v-i.  we  use  the 
notation  V(t’)  and  V''(^,u)  whenever  it  is  necessary. 

Assuming  that  all  the  entries  of  V{z)  are  rational  functions  of  r  vanishing  at  infin¬ 
ity  with  poles  having  negative  real  part,  and  taking  V  as  an  operator  on 
(1  <  p  <  oo),  the  spectral  radius  p{V)  =  linifc^oo  satisfies  ([6]) 

(12)  p{V)=mAxp{V(z)). 

H€Z>0 

Note  that  .S'(0)  and  F(0)  are  respectively,  the  smoothing  and  multigrid  two-grid  cycle 
operators  for  the  corresponding  steady-state  problem  (or  static  problem)  LhU  =  /. 


Theorem  1.  The  spectral  radius  for  the  waveform  two-grid  V-cycle  operator 
V  =  V(i;)  for  equation  (5)  satisfies 

(13)  lim  |p(V(u))-p(l/(0,u))l  =  0 

v—*oo 

for  either  daiaped  Jacobi  or  red-black  Gauss-Seidel  smoothings. 

Proof.  We  shall  prove  the  case  of  the  red-black  Gauss-Seidel  relaxation  onl}'.  The 
proof  for  the  damped  Jacobi  relaxation  follows  similarly.  Let  CG{z)  denote  the  Laplace 
transform  of  coarse  grid  correction  operator, 

CGiz)  =  /-/^(2-t-LH)-V"(r-f  U) 

=  C6'(0) -k  GG\(z) 
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(14) 


with 


(15)  CG^(z)  =  +  LH)-'(L-„'lj:U  -  /»"). 

Since  the  spectral  radius 

p{S'’^{z)CG{z)S'^'{z))  =  p{CG{z)S^{z)),  v  =  v,  +  1-2,  for  all  c, 
for  simplicity,  we  rewrite  the  two-grid  operator  (II)  as 


(16) 


Viz,v)  =  CG{z)S'-iz) 

=  CG{0)S'’{z)  +  CG^{z)S^iz). 


As  noted  in  Section  2,  the  red-black  Gauss-Seidel  smoother  (8) 


(17) 


•5”(^)  =  (its) 


2v 


o 

-- 

(^)  BC'’-'  ' 

Lo 

C'’ 

C  =  RB, 


,  2v 


v2v 

’  0 

) 

0 

0 

Let  2  :=  z/d,  we  have 


^"("7^)  =  7rTTT5;;mt^)  +  A(2,n) 


with 

(18) 


(i  +  I)'" 

A(2,u)  =  F(‘'>(^)A(2,i;), 


=  - 


(z-f  1)2>'’ 


A(2,i;)  =  CG(0) 


0 

0  0 


0  {\  +  ~)Br 

0  C" 


Given  f  >  0,  there  is  a  matrix  norm  1|  •  [j,,  such  that  (see  [4],  pp.297) 


p{V{z,v))  < 


'{z+l) 


2v 


V(0,u)||.-f  ||A(2,e)l|, 


and 


(2  +  1) 


2v 


V^(0,n))  +  e. 


Therefore 


1 


p(V(u))  <  max  '  p{V{Q,  v))  +  c  +  ma^xJ|A(z,  r)ll,. 
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which  leads  to 


(19)  0  <  /j(V(u))  -  <  f  +  max  ||A(c,v)||<. 

He2>0 

The  matrix  function  A(5,i;)  is  analytic  in  Rez  >  0,  continuous  on  Rez  =  0,  and 
A(oo,u)  is  bounded  for  all  u  >  1  (because  .S'''(0)  — >  0  as  i;  — >  od  [7]),  so  does  its  e-norm 
||A(2,t;)||,,  thus 

(20)  max  |1A(2,u)||,  <  c(e), 

Hez>Q 

where  c(f)  is  a  function  of  e  independent  of  z  and  v.  Combining  Eq.(18)  and  Eq.(20), 
we  iiave 

max  ||A(.,.)||,  <  maxJFH,,), .  .  „„ 

— >  |F^''^(0)|  •  c(e)  =  0,  as  u  — +  oo.  □ 

4.  Fourier-Laplace  analysis.  Fourier  analysis  has  been  used  to  provide  exact 
convergence  rate  of  multigrid  iteration  for  some  steady-state  problems  [1].  In  order  to 
analyze  the  convergence  rate  of  waveform  multigrid  iteration,  first,  the  Laplace  trans¬ 
form  is  introduced  in  time  variable  to  convert  a  m-dimensional  semi-discrete  linear 
differential  equation  with  time-independent  coefficients 

(21)  -t-  LfjU  =  / 

u(0)  =  0 

to  equivalent  algebraic  equations  defined  over  the  right  half  complex  plane 

(22)  (2  Lh)u{z)  =  /(z),  Rez  >  0. 

Then  multigrid  method  applied  in  space  only  is  analyzed  by  considering  its  action  on 
each  of  these  equations.  This  is  done  using  Fourier  analysis  for  each  equation  in  (22), 
and  combining  the  results  to  obtain  the  convergence  rate  of  waveform  multigrid  iteration 
for  Eq.(21).  In  this  section  a  Fourier  Laplace  analysis  is  applied  to  the  waveform  two 
grid  V-cycle  iteration  described  in  Section  3. 

Let  us  begin  with  the  error  formula  (10) 

g(.)  3, 

Its  Laplace  transform  is 

(23)  e<’)(2)  =  l/(2)e<’-’)(2),  Rez  >  0. 
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Let 


9  =  =  0’(mod  tt),  j  =  2,---r\ 

and  X{9),  V{0,z)  denote  the  Fourier  modes  on  grid  h  and  the  symbol  of  V(z),  i.e., 

X{0)  =  \exp{iB^x ! h),  •  •  - ,  exp{i6^ 

V{9,z)  =  -)]2'”x2’” 

satisfying 

(24)  V{z)X{e)  =  X{d)Vie,z)  forall^. 

Then 

(25)  piV{2))  =  supp{V{e,z))  =  m&x  p{V {9,  z)). 

e  o<l®M<f 

The  symbol  matrix  V{9,  z)  is  of  order  2”*,  a  tiny  matrix  comparing  to  V{z).  Its  spectral 
radius  can  be  calculated  accurately  for  each  given  9  and  z.  Therefore  p(V),  the  conver¬ 
gence  rate  of  the  waveform  two  grid  iteration,  is  obtained  by  computing  p{V[9,z)  over 
0  <  <  f  and  Rez  =  0: 

/?(V)  =  max  max  /)(F(0,2)). 

''  flej=0  0<|e'|<f 

Example.  Consider  the  heat  equation  on  a  square  fl  =  (0,  tt)  x  (0,  tt)  with  Dirichlet 
boundary  conditions 

Ui  —  Au  =  /,  (<,x)  €  (0,  oo)  X  n 

u  =  g,  {t,x)  G  [0,oo)  X  dfl,  u(0,x)  =  Uo(x),  x  E  fl. 

Let  L/i  correspond  to  the  five  point  Laplacian 


The  Laplace  transform  of  the  two  grid  iteration  operator  V  is  given  by  (see  Eq.(16)) 
(27)  V(z)  =  CC(z)S^(z) 

CG(z)  =  /  -  /,)(z  +  Lh)-'  lUiz  T  U) 
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with 


1 


and  the  corresponding  smoother  ^{2).  Since  the  space  spanned  by  the  Fourier  mode 
X(0)  is  invariant  under  each  of  operators  in  Eq.(27),  we  have 

(28)  V{e,z)  =  {I  -  i^„{6){z  +  L„{26))-^ll^{e)(z+  h{0)))S^e,z),  Rez  >  0. 


where  ~’s  indicate  the  matrix  symbols  and  H  =  2h  is  assumed.  For  this  particular 
example,  the  fine  grid  operator  is  represented  as 


1,(0)  = 


L,{0^) 


U{0^)  J 


with 


0  =  {0U02),  Ui0)  = 

the  coarse  grid  operator  as 


4(szn^(^,/2)  +  sm^(^,/2)) 


LHi20)  =  {LH{20^)]. 

The  bi-linear  interpolation  is  chosen  for  and  the  restriction  operator 
Their  matrix  symbols  are 


Smoothers  for  the  damped  Jacobi  and  the  red-black  Causs-Seidel  relaxation  are  repre¬ 
sented  as  (see  Eq.(7)-(8)) 


5.y(0,c)  = 


S{0\z) 


Si0\z]  J 


—l~{d-u;U{e)),  0<u;<l; 

zu  +  a 


Sgs[0.z)  = 


5',.(^,z)  0 

0  Sk{0,z) 


o  ^  _  M  ^  -(1  + 


1  -f-  «  — ( 1  -F  fl) 

1  —  a  —(1  —  a) 
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Table  1 

Damped  Jacobi  Waveform  Muliigrt  =  2/3) 


- 1 

f  1 

1 

2 

3 

4 

5 

6 

7 

8  9  10 

p{V(0) 

.6626 

.4390 

.2909 

.1927 

.I486 

.1269 

,1124 

.1004  .0908  .0838 

p(V'(0,r)) 

,6626 

.4390 

.2909 

.1927 

13.32 

.11.52 

,1003 

.0881  .0793  ,0728 

Table  2 

Red-Black  Gauss-Setdel  Waveform  Mullignd 


1 

2 

3 

4 

5 

6 

( 

8 

9 

10 

p(V(i.)) 

.2439 

.1.567 

.1132 

.0861 

.0737 

.0614 

.0.548 

.0473 

.''too 

.0390 

/ifVTO.i;)) 

.24,39 

,1517 

1046 

.0765 

,0669 

.0,543 

.0474 

,0404 

.0.368 

.0336 

“  =  2(1  ‘  =  20^' 

Tables  1  and  2  present  computed  p(V(i>)),  the  spectral  radius  of  waveform  multi¬ 
grid  operator,  and  /j(V^(0,n))  of  corresp  tiding  steady-state  multigrid  operator  for  the 
damped  Jacobi  and  the  red-black  Gauss-Seidel  relaxation.  As  it  shows,  the  estimated 
convergence  rates  of  waveform  multigrid  iteration  are  virtually  as  good  as  those  for 
steady-state  problems.  Since  extensive  research  has  been  done  in  mulligrid  methods 
for  steady-state  problems  [1]  [3]  [10],  this  analytic  comparison  is  very  useful  in  pre¬ 
dicting  the  performance  of  waveform  multigrid  methods.  Although  our  discussion  was 
done  only  for  the  model  heat  equation,  one  may  expect  t  e  same  performance  of  the 
waveform  multigrid  iteration  as  that  of  sieady-state  one  in  wide  cicisses  of  applications. 
The  treatment  of  general  problems,  e.g.,  non-constant  coefficients  or  general  domains,  is 
done  in  a  similar  way  using  frozen  coefficients  argument  which  is  applicable  for  smooth 
coefficient  problems.  Its  theoretical  rigorous  justification  is  involved  and  needs  the  use 
of  pseudo-difference  calculus. 

5.  Comparison  and  Conclusion.  We  have  shown  that  the  estimated  conver 
gence  rate  of  waveform  multigrid  iteration  is  almost  undistinquishable  from  that  of 
steady-state  one.  Now,  we  compare  the  experimental  results  of  the  convergence  rate  to 
the  analytic  ones  discussed  in  last  section. 

Consider  the  two-dimensional  heat  equ.ition  on  il  =  (O.tt)  x  (0,  tt) 

U(  —  All  =  /,  (c  x)  €  (0,  /  /]  X  n, 

nit,x)  =  0.  (Cx)  e  [0,  </]  X  dn, 

ii(o,  j)  =  uo(x).  X  t  n, 

where  uq  was  generated  randomly  to  excite  all  possible  Fourier  modes.  The  space 
derivatives  were  discretized  by  central  differences  with  uniform  fine  grid  size  hj..  For 
time  integration,  as  in  [7],  trapezoidal  rule  wa,s  used.  7’he  step  size  was  always  chosen 


Table  3 

Damped  Jacobt  Waveform  Multigrid 


analytic  p{V(t;)) 

measured  p(V(j.')) 

worst  case 

avareage  case 

1 

.6626 

.66.57 

.6403 

2 

.4390 

.4401 

.4160 

3 

.2909 

2964 

.2693 

4 

.1927 

.1890 

1836 

Table  4 

Red-Black  Gauss-Setdel  Waveform  Multigrid 


V 

analytic  p{V{v)) 

measured  p(V(v)) 

iy/r;n(2v  -  1 )  -1-  O(lil) 

worst  case 

avareage  Ccise 

I 

.2439 

.2.560 

.2282 

.■2b00-i-O(hi) 

2 

.1567 

.1137 

.0740 

.mb+0{til) 

3 

.1132 

.0761 

.0524 

.1295-(-0(/i2) 

4 

.0861 

.0670 

.0437 

.1110+0(^2) 

Bs  ht  =  .01  on  [0,</].  Note,  the  efficiency  of  time  integration  is  not  the  concern  of  this 
paper. 

The  experiments  were  run  for  two-grid  V-cycle  with  v  damped  .Jacobi  (u)  =  2/3) 
or  red-black  Gauss-Seidel  relaxations.  The  fine  grid  mesh  size  in  space  was  =  7r/n, 
n  =  8,16,32,64.  Although  the  spectral  radius  p{V{v))  was  studied  on  entire  time 
interval  t  G  [0,oo),  finite  interval  [0,tf]  had  to  be  used  in  experiments.  For  each  set 
of  tests,  we  used  both  tf  =  [  and  tf  =  10.  However,  we  found  that  the  measured 
convergence  rates  did  not  depend  on  the  size  of  time  interval. 

To  measure  the  spectral  radius  p(V(w)),  we  used  the  asymptotic  ratio  of  the  defects 
(see  Eq.(9)) 


max(  ||4"^’^|i2 
max,  ||4'^1|2 

Using  the  matrix  split  of  Lh,  the  calculation  of  the  derivatives  in  the  defects  can  be 
avoided.  Pecause  of  extensive  computations  involved,  ||d|,'*|l2  was  evaluated  only  at 
t  =  tf  since  the  convergence  of  the  waveform  iteration  is  determined  by  the  error  at 
t  =  tf. 

Tables  3  and  4  present  the  comparison  results.  They  show  that  the  analytic  con¬ 
vergence  rates  obtained  via  Fourier- Laplace  analysis  almost  coincide  with  the  real  ones, 
a  result  evidenced  for  the  steady-state  problems  [1].  As  a  reference,  Table  4  also  lists 
values  of  |^7/o(2u  —  1),  the  bounds  of  p{V{v))  proved  for  one-dimensional  problem  by 
Lubich  and  Ostermann  [6].  They  conjectured  the  bounds  for  that  of  two-dimensional 
heat  equation  as  ~y^Tjo{2v  —  1)  -f-  0{hl).  Observing  our  analysis  in  Section  4,  with 
standard  discretization  of  L,  nine-point  restriction  and  bi-linear  interpolation,  /)(V(r)) 
is  independent  of  the  grid  size  hx-  The  major  difference  of  our  approach  to  the  one 
in  [6]  is  that,  instead  of  using  eigenvectors  of  Lh,  we  used  Fourier  modes,  which  are 
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much  easier  to  manipulate,  extendible  to  high-dimensional  problems  and  wide  class  of 
applications.  Most  important  of  all,  this  approach  is  able  to  give  the  exact  convergence 
rates  for  special  model  problems  and  sharp  estimates  for  general  problems. 
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